/*
Resampling method encapsuated in a class. 
Input: STL vector of type double
Output: STL vector of integer indicating resampling indecies
*/
#pragma once
#include <time.h>
#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <vector>
#include "roundup.h"
#include <math.h>
#include <algorithm>
#define DEMO_PRINT 0
#define PROGRESS_PRINT 1
using namespace std;


class Resampling {
public:
	Resampling(void);
	~Resampling(void);
	vector<int> resample(vector<double> weights) const;
	vector<int> resample_monte_carlo(vector<double> weights) const;
	void resample_naive_test(int length = 3)const;
	double compute_correlation(vector<double> vector1, vector<double> vector2) const;
private:
	int region(double particle,vector<double> weights) const;
protected:

};

Resampling::Resampling(){
}

Resampling::~Resampling(){
}

int Resampling::region(double particle,vector<double> weights) const{
	int index = 0;
	double sum = 0;
	for (int i = 0; i < weights.size(); i++){
		sum += weights[i];
		if (sum>particle){
			index = i;
			break;
		}
	}
	return index;
}

vector<int> Resampling::resample_monte_carlo(vector <double> weights) const{
	
	vector<double> new_sample_number;
	vector<int> new_sample_indeces;
	int index_to_push, rand_buf;
	double particle;

	if(DEMO_PRINT) cout<<" Input------Rand-------Particle-------Resample"<<endl;
	for (int i = 0; i < weights.size(); i++)
	{
		if(DEMO_PRINT) cout<<i<<"	"<<weights[i]<<"	";
		//RAMDONLY GENERATE SAMPLE
//		srand(time(NULL));
		rand_buf = rand();
		particle = (double)(rand_buf % 10000 + 1.0)/10000.0;	//0-1		//EVEN IF rand() IS UNIFORMAL, THIS MIGHT NOT BE UNIFORM DISTRIBUTED
		if(DEMO_PRINT) cout<<particle<<"	";
		//COMPUTE WHICH AREA IT FELL INTO 
		index_to_push = region(particle, weights);		
		if(DEMO_PRINT)  cout<<" "<<index_to_push<<endl;
		//REPORT INDEX
		new_sample_indeces.push_back(index_to_push);
	}
	if(DEMO_PRINT)  cout<<endl;
	sort(new_sample_indeces.begin(), new_sample_indeces.end());
	return new_sample_indeces;
}

void Resampling::resample_naive_test(int length)const {
//	int length = 3;
	vector<double> test_vector;
	vector<int> distribution;
	vector<double> distribution_normalized;

	for (int i = 0; i < length; i++)
	{
		test_vector.push_back(rand());
		distribution.push_back(0);
	}
	if(PROGRESS_PRINT) cout<<"INPUT FINISHED"<<endl;

	//NORMALIZE RANDOMLY GENERATED PROBABILTY
	int sum = 0;
	for (int i = 0; i < test_vector.size(); i++)
	{
		sum += test_vector[i];
	}

	for (int i = 0; i < test_vector.size(); i++)
	{
		test_vector[i] = test_vector[i]/sum;
	}
	if(PROGRESS_PRINT) cout<<"NORMALIZATION FINISHED"<<endl;


	sum = 0;
	vector<int> resample_indecies = resample_monte_carlo(test_vector);
	for (int i = 0; i < resample_indecies.size(); i++)
	{
		if(DEMO_PRINT) cout<<" "<<resample_indecies[i];
		distribution[resample_indecies[i]] += 1;
	}
	if(DEMO_PRINT) cout<<endl<<"Distribution:"<<endl;
	if(PROGRESS_PRINT) cout<<"MONTE_CARLO SIMULATION FINISHED"<<endl;

	double correlation = 0, euclidean_distance = 0;
	for (int i = 0; i < distribution.size(); i++)
	{	
//		if(PROGRESS_PRINT) cout<<"Computing Euclidean Distance for "<<i<<endl;
		euclidean_distance += pow(((double)distribution[i]/distribution.size() - (double)test_vector[i]),2.0);
		distribution_normalized.push_back((double)distribution[i]/distribution.size());
		if(DEMO_PRINT) cout<<i<<"		"<<test_vector[i]<<"		"<<distribution_normalized[i]<<endl;
	}
	if(PROGRESS_PRINT) cout<<"Euclidean_distance = "<<euclidean_distance<<"	";
	
	correlation = compute_correlation(test_vector,distribution_normalized);
	if(PROGRESS_PRINT) cout<<"Correlation = "<<correlation<<endl;

	cin.get();

}

double Resampling::compute_correlation(vector<double> vector1, vector<double> vector2) const{
	double mean1 = .0;
	double mean2 = .0;
	double s1, s2, cov, correlation;
	double buffer1 = .0;
	double buffer2 = .0;
	double buffer0 = .0;
	//COMPUTE MEAN
	for (int i = 0; i < vector1.size(); i++){		//THEY HAVE TO HAVE SAME SIZE
		buffer1 += vector1[i];
		buffer2 += vector2[i];
	}
	mean1 = buffer1/vector1.size();
	mean2 = buffer2/vector2.size();

	//COMPUTE STANDART_DEVIATION
	buffer1 = .0; buffer2 = .0;
	for (int i = 0; i < vector1.size(); i++){		//THEY HAVE TO HAVE SAME SIZE
//		if(PROGRESS_PRINT) cout<<"Computing correlation for "<<i<<endl;
		buffer1 += pow((vector1[i] - mean1),2.0);
		buffer2 += pow((vector2[i] - mean2),2.0);
		buffer0 += (vector1[i] - mean1)*(vector2[i] - mean2);
	}
	s1 = sqrt((buffer1)/(vector1.size() - 1));
	s2 = sqrt((buffer2)/(vector2.size() - 1));
	cov = buffer0/(vector1.size() - 1);

	correlation = cov/(s1*s2);
	return correlation;
}
